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1  Summary 


The  present  project  has  been  initiated  with  the  aim  to  produce  a  solution  methodology  for  high- 
order  discretization  of  compressible  flow  problems  that  substantially  meets  industry  standards  in 
terms  of  efficiency  and  robustness.  The  approach  is  two-pronged:  Both  discretization  methods  and 
solution  techniques  for  the  arising  nonlinear  systems  of  equations  are  addressed. 

This  report  focuses  on  the  results  achieved  during  the  most  recent  no-cost  extension  of  the 
project.  It  is  mostly  centered  around  novel  discretization  methods,  in  particular  hybridized  DG 
schemes.  During  the  current  funding  phase,  we  have  further  developed  our  hybridized  DG  dis¬ 
cretization  method  for  compressible  flow  simulation. 

Hybridization  has  been  identified  a  potential  breakthrough  technology  allowing  one  to  reduce 
the  number  of  globally  coupled  degrees  of  freedom  (DOFs)  very  significantly:  Using  polynomials 
or  order  m  with  conventional  discretization,  globally  coupled  DOFs  grow  as  and  m?  in  two  and 
three  dimensions,  respectively,  while  for  a  hybridized  discretization  this  reduces  to  and  m.  This 
is  obviously  very  significant  for  both  computational  efficiency,  and  storage  requirements,  which  has 
been  recognized  as  a  major  bottleneck  for  implicit  solution  methodologies  with  high-order  methods. 
Among  the  highlights  that  will  be  exposed  in  more  detailed  in  the  technical  section  of  this  report 
are 


•  Implementation  and  validation  of  shock  capturing  capability 

•  Target-based  hp-adaptation  techniques 

During  the  most  recent  funding  period,  research  partially  funded  by  this  project  has  been 
presented  at 

•  World  Congress  on  Computational  mechanics,  July  2012,  Sao  Paulo,  Brazil 

•  ECCOMAS  congress  on  computational  methods  in  applied  sciences  and  engineering,  Septem¬ 
ber  2012,  Vienna,  Austria 

•  International  Workshop  on  High-Order  CFD  Methods,  May  2013,  Cologne,  Germany 

An  updated  list  of  published  results,  covering  the  entire  funding  period,  which  acknowledge  support 
from  the  current  grant  is  given  below: 

•  hybridized  DG  schemes,  including  target-based  adaptation  via  adjoint  equation  [6,  7,  8] 

•  Stable  high-order  Spectral  Difference  method  for  hyperbolic  conservation  laws  on  triangules  [1] 

•  Multilevel  methods  for  solution  of  the  Euler  equations  [5]; 

•  Matrix-free  solution  methods  for  implicit  relaxation  schemes  in  response  to  the  extreme  stor¬ 
age  requirement  of  implicit  relaxation  methods  for  high-order  discretization  [5,  3,  2]; 

•  Equivalence  between  Spectral  Difference  (SD)  and  nodal  Discontinuous  Galerkin  (DG)  schemes  [4] 
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2  Governing  Equations 

We  consider  systems  of  partial  differential  equations 

V  ■  (iciw) -fv{w,'Vw))  =  s{w,'Vw)  (1) 

with  convective  fluxes  fc  :  M™'  — )■  and  diffusive  fluxes  ft,  ;  x  — )■  and  a  state- 

dependent  source  term  s  :  x  — )•  M™.  Potentially,  some  of  these  quantities  could  be  zero. 

We  denote  the  spatial  dimension  by  d  and  the  number  of  conservative  variables  by  m. 


2.1  Two-Dimensional  Euler  Equations 

The  Euler  equations  are  given  in  conservative  form  as 

V  •  f,H  =  0  (2) 


with  the  vector  of  conserved  variables 


w  =  {p,pv,Ey 


(3) 


where  p  is  the  density,  v  is  the  velocity  vector  v  :=  {vx,Vy)'^ ,  and  E  the  total  energy.  The  convective 
flux  is  given  by 


fc  =  {pv,pld  -h  V  (g)  v,v{E  +  p))'^  . 

Pressure  is  related  to  the  conservative  flow  variables  w  by  the  equation  of  state 

p={j-l)(^E-  ^pv  •  V 

where  7  =  Cp/c„  is  the  ratio  of  specific  heats,  generally  taken  as  1.4  for  air. 
Along  wall  boundaries  we  apply  the  slip  boundary  condition 

Vn{w)  :=  V  •  n  =  0. 

We  also  define  a  boundary  function  which  satisfies  Vn{wdn{w))  =  0  as 

In  1  —  — n  n 

wdn{w)  = 


0 


vj  -L  — 

0  —UxUy  ^  —  ny  0 


w. 


(4) 


(5) 


(6) 


(7) 


\0  0  0  1/ 

At  the  far-held  can  be  realized  with  the  aid  of  characteristic  upwinding  (exposition  omitted) 


2.2  Two-Dimensional  Navier-Stokes  Equations 

The  Navier-Stokes  equations  in  conservative  form  are  given  by 

V  •  (fc(tc)  -  fc(te,  Vrc))  =  0.  (8) 

The  convective  part  fc  of  the  Navier-Stokes  equations  coincides  with  the  Euler  equations.  The 
viscous  hux  is  given  by 

fy  =  (0,  r,  TV -I- A:Vr)^  .  (9) 
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The  temperature  is  defined  via  the  ideal  gas  law 


T  = 


k-Pi  \p  2  J  (7  -  l)c„  p 


(10) 


where  Pr  =  ^  is  the  Prandtl  number,  which  for  air  at  moderate  conditions  can  be  taken  as  a 
constant  with  a  value  of  Pr  ~  0.72.  k  denotes  the  thermal  conductivity  coefficient.  For  a  Newtonian 
fluid,  the  stress  tensor  is  defined  as 


T  =  Vv  +  (Vv)^  -  ^  (V  •  v)  Id  )  . 

O 


(11) 


The  variation  of  the  molecular  viscosity  //  as  a  function  of  temperature  is  determined  by  Suther¬ 
land’s  law  as 

-75 

with  Cl  =  1.458  X  10“®  and  C2  =  110.4K. 

msv  K 

Along  wall  boundaries,  we  apply  the  no-slip  boundary  condition,  i.e. 


V  =  0 


(13) 


with  corresponding  boundary  function 


waniw)  =  {p,0,0,Ey 


(14) 


Furthermore,  one  has  to  give  boundary  conditions  for  the  temperature.  In  the  present  work  we  use 
the  adiabatic  wall  condition,  i.e. 

Vr-n  =  0  (15) 


Combining  both  no-slip  and 
flux,  namely 


adiabatic  wall  boundary  conditions,  gives  a  condition  for  the  viscous 


fv,dn  {wdn,<idn) 


0  Til  T21  0\^ 

0  Ti2  T22  oj 


(16) 


3  Discretization 

3.1  Notation 

We  tesselate  the  domain  into  a  collection  of  non-overlapping  elements,  denoted  by  Th,  such  that 
Uii'eT^  K  =  Pi.  For  the  element  edges  we  consider  two  different  kinds  of  sets,  dTh  and  F/j,  which 
are  element-oriented  and  edge-oriented,  respectively. 

dTh  :=  {  dK\dPl  :  KeTh},  (17) 

Fh  :=  {  e  :  e  =  K  D  K' for  K,  K' G  Th',^easd-iie)  ^  0  }.  (18) 

The  first  is  the  collection  of  all  element  boundaries,  which  means  that  every  edge  appears  twice. 
The  latter,  however,  includes  every  edge  just  once.  The  reason  for  this  distinction  will  become  clear 
later.  Please  note  that  neither  of  these  sets  shall  include  edges  lying  on  the  domain  boundary;  the 
set  of  boundary  edges  is  denoted  by  r)(. 
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We  denote  by  11^ (D)  the  set  of  polynomials  of  degree  at  most  p  on  some  domain  D.  We  will 
need  discontinuous  function  spaces  for  the  domain  and  the  mesh  skeleton; 


Yh  =  {v£  (D) 

v\k 

£  np^(i^). 

K  £  Th}^^^ 

(19) 

Wh  =  {w£  (D) 

w\x 

G  np^(i^). 

K  £  nr 

(20) 

Mh  =  {p£  (F;,) 

P\e 

G  TT>^{e), 

e  GF4-. 

(21) 

Thus,  V  £  Vh,  w  £  Wh  and  p  £  Mh  are  piecewise  polynomials  of  degree  p  which  can  be  discontin¬ 
uous  across  edges  (for  v,  w)  or  vertices  (for  p),  respectively. 

Usually,  the  polynomial  degree  between  elements  and  interfaces  does  not  vary.  In  the  case 
of  varying  polynomial  degrees  px-  and  Px+i  we  choose  the  polynomial  degree  for  the  interface 
e  =  K~  n  K'^  as  pe  =  max  {px-  ■,  Pk+  }  ■ 

We  distinguish  between  element-oriented  inner  products  (defined  with  Th)  and  edge-oriented 
inner  products  (defined  with  Th) 


{v,w)-^^:=  /  uu;dx, 

K&Th 

(v,w).^^  ■■=  Y  /  V- wdx. 

(22) 

K&Th 

{v,w)y^  ■=Y  '^wda. 

eer^ 

(23) 

3.2  Weak  Formulation 

We  can  rewrite  general  convection-diffusion  equations  as  a  first-order  system  by  introducing  an 
additional  unknown  representing  the  gradient  of  the  solution 

q  =  Vw  ^  ^ 

(24 

V  •  (fc  (w)  -  f„  {w,  q))  =  s  {w,  q) . 

By  multiplying  the  strong,  mixed  form  (24)  with  appropriate  test  functions  {Th,  Ph)  £  h  x  Wft 
and  integrating  by  parts,  we  obtain  a  standard  DG  discretization  in  mixed  formulation  of  the 
problem,  i.e.; 

Find  {<\h,Wh)  £VhxWh  s.t.  \/{Th,Ph)  £^hXWh 


0=  {(\h,Wh-,Th,'Ph)  (25) 

:=  (t/i,  cih)r,,  +  (^  •  Wh)rH  ~  ■  n,  w)q^^  (26) 

-  {Viph,T{wh)  -  iv{wh,c^h))r,.  -  {^h,s{wh,cih))x  +  \  Th,!c  -  fv)  (27) 

{cih,Wh-,Th,Ph)  +  ^fh%  {cih,Wh-,Ph)  ■  (28) 


Here  the  numerical  trace  w  and  the  numerical  fluxes  fc,  fv  have  to  be  chosen  appropriately  to 
define  a  stable  and  consistent  method.  Furthermore,  the  boundary  conditions,  here  denoted  by 
■^hdn  (H/o  Wh',Th,  Ph),  have  to  be  discretized  appropriately. 

In  contrast  to  a  DG  discretization,  where  the  numerical  trace  w  is  defined  explicitly  in  terms  of 
Wh  and  c[h,  it  is  treated  as  an  additional  unknown  in  an  HDG  method.  This  additional  unknown 
is  called  Xh  and  has  support  on  the  skeleton  of  the  mesh  only  In  order  to  close  the  system  the 
continuity  of  the  numerical  fluxes  across  edges  is  required  in  a  weak  sense,  resulting  in  a  third 
equation. 
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The  weak  formulation  of  the  hybrid  system,  comprised  of  equations  for  the  gradient  q/j,  the 


solution  itself  Wh  and  its  trace  on  the  mesh  skeleton  is  then  given  by: 

Find  (q;,,  Wh,  Xh)  G  ;=  x  IT/^  x  Mh  s.t.  V(rfe,  (ph,  fih)  G 

0=  JVh{qh,Wh,Xh;Th,(ph,Ph)  (29) 

:=  (t'/i)  ^  ~  dTh  (^1^) 

-  {Vph,^c{wh)  -  f„(w/i,q/i))-T-  -  {ph,s{wh,cih))r^  +  {ph,fc  -  fv)  (31) 

'*  "'  \  /  dTh 

+  Xfh,dQ  (q/i,  Wh;Th,  Ph)  +  A6i.sc  (q/i,  Wh;  Ph)  (32) 

+  Uh,  \fc  -  7J  )  .  (33) 


The  terms  tested  against  Th  and  ph  are  called  local  solvers,  meaning  they  do  not  depend  on  the 
solution  within  neighboring  elements  but  only  on  the  trace  of  the  solution  which  is  approximated  by 
Xh-  The  coupling  between  elements  is  then  introduced  by  weakly  enforcing  the  normal  continuity 
of  the  numerical  fluxes  across  interfaces. 

We  choose  numerical  fluxes  comparable  to  the  Lax-Friedrich  flux  and  to  the  LDG  flux  for  the 
convective  and  diffusive  flux,  respectively,  i.e. 

fc  (Aft,  Wh)  =  fc  (Aft)  •  n  -  ac  (Aft  -  Wh)  (34) 

fy  {Xh,  Wh,  qft)  =  iy  {Xh,  qft)  •  n  +  ay  {Xh  -  Wh)  (35) 

which  can  be  combined  into 

fc-fv  =  (fc  (Aft)  -  ft,  (Aft,  qft))  •  n  -  (ac  +  a„)  (Aft  -  Wh)  -  (36) 

The  stabilization  introduced  can  be  given  by  a  tensor;  in  our  work,  however,  we  restrict  ourselves 
to  a  constant  scalar  a  =  ay  +  ay  which  seems  to  be  sufficient  for  a  wide  range  of  test  cases. 

3.2.1  Boundary  Conditions 

In  order  to  retrieve  an  adjoint-consistent  scheme,  special  care  has  to  be  taken  when  discretizing 
the  boundary  conditions  (see  [7]).  The  boundary  conditions  have  to  be  incorporated  by  using  the 
boundary  states  wqq  {wh)  and  gradients  qgo  {wh,  qft),  i.e. 

A/ft, an  (qft,  Wh;  Th,  Ph)  ■=  {rh  •  n,  Wdn)Yb 

h 

+  {Ph,  (fc  {wan)  -  fc  {wan,  qao))  •  n)p6  . 

We  would  like  to  emphasize  that  Aft  does  not  occur  in  this  boundary  term. 

3.2.2  Shock-Capturing 

We  adopt  a  shock-capturing  term,  where  an  artificial  viscosity  term,  given  by  V  •  (e  {w,  Vrc)  Vtc), 
is  used.  The  viscosity  e  is  given  by  the  Ll-norm  of  the  strong  residual  V  •  fc(rc)  in  every  element. 
In  order  to  accelerate  the  convergence  of  this  term  to  zero  with  mesh  refinement,  it  is  premultiplied 
with  an  effective  mesh  size  Jik  '■=  —■  The  latter  resembles  the  actual  resolution  within  an  element. 

PK 

Furthermore,  a  user-defined  factor  cq  is  introduced  which  can  be  reliably  tuned  for  a  rather  large 
range  of  test  cases.  Finally,  the  artificial  viscosity  is  given  by 
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where  the  strong  residual  is  given  by 

m 

d(y;)  :=  ^|(V-fc(u;))i|  .  (38) 

i=\ 

In  the  discretization  of  this  shock-capturing  term  the  interface  integral  is  neglected  so  that  only 
the  volume  contribution  is  considered,  i.e. 

A4,sc  [Wh]  ^h)  ■=  e  {Wh,  Vwh)  ^Wh)-Y^  ■  (39) 

This  obviates  the  need  for  introducing  q/j  in  purely  convective  problems  (e.g.  the  compressible 
Euler  equations).  In  the  viscous  case,  where  the  gradient  is  explicitly  given,  can  be  replaced 
by  q/j  yielding 

A/'/i, sc  (q^,  Wh]  (Ph)  ■=  (Vv9ft,  e  {wh,  qh)  q/i)r^  •  (40) 

Please  note,  that  this  term  enters  only  the  local  part  of  the  discretization. 

3.3  Relaxation 

We  solve  the  nonlinear  system  of  equations  that  defines  the  HDG  method,  with  a  damped  Newton- 
Raphson  method.  An  artihcial  time  is  introduced,  and  we  solve  at  each  iteration  index  n, 

(^Ph,  +-^h  K]  =  -^h  «;yft)  Vyfe  G  Xft.  (41) 

Please  note  that  by  choosing  At”  — )■  oo,  a  pure  Newton-Raphson  method  is  obtained.  Usually  the 
time  step  is  kept  finite  for  a  few  initial  steps  to  ensure  stability.  As  soon  as  the  residual  is  lower 
than  a  certain  threshold,  i.e.  the  current  approximation  x))  is  thought  to  be  sufficiently  close  to 
the  solution  x/i,  we  let  the  time  step  go  towards  infinity. 


3.4  Hybridization 


Using  an  appropriate  polynomial  expansion  for  hq/i,  5wh  and  5\h,  the  linearized  global  system  is 
given  in  matrix  form  as 


'  A  B  R' 

■  SQ  ■ 

■  F  ■ 

CDS 

6W 

= 

G 

L  M  N 

SA 

H 

(42) 


where  the  vector  [5Q^  dlU,  dA]^  contains  the  expansion  coefficients  of  hx/j  with  respect  to  the  chosen 
basis. 

In  order  to  carry  on  with  the  derivation  of  the  hybridized  method,  we  want  to  formulate  that 
system  in  terms  of  dA  only.  Therefore  we  split  it  into 


■  A 

B  ' 

■  5Q  ' 

■  F  ' 

■  R  ' 

C 

D 

SW  _ 

G 

S 

<5A 


and 


■  L  M  ■ 


5Q 

5W 


N5K  =  H. 


Substituting  Eq.  (43)  into  Eq.  (44)  yields  the  hybridized  system 


N  -[  L  M  ' 


■  A 

B  ' 

-1 

■  R  ' 

C 

D 

S  _ 

5A  =  H  -  [  L  M  ' 


■  A 

B  ' 

-1 

■  F  ' 

C 

D 

G 

(43) 

(44) 


(45) 
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The  workflow  is  as  follows;  First,  the  hybridized  system  is  assembled  and  then  being  solved  for 
5 A.  Subsequently,  5Q  and  SW  can  be  reconstructed  inside  the  elements  via  Eq.  (43).  It  is  very 
important  to  note  that  it  is  not  necessary  to  solve  the  large  system  given  by  Eq.  (43).  In  fact,  the 
matrix  in  Eq.  (43)  can  be  reordered  to  be  block  diagonal.  Each  of  these  blocks  is  associated  to 
one  element.  Thus,  both  the  assembly  of  the  hybridized  matrix  in  Eq.  (45)  and  the  reconstruction 
of  SQ  and  5W  can  be  done  in  an  element-wise  fashion.  In  order  to  save  computational  time,  the 
solutions  to  Eq.  (43)  can  be  saved  after  the  assembly  of  the  hybridized  system  and  reused  during 
the  reconstruction  of  5Q  and  5W . 

The  hybridized  matrix  is  a  nj  x  nj  block  matrix,  where  nj  =  |r/i|  is  the  number  of  interior 
edges.  In  each  block  row  there  is  one  block  on  the  diagonal  and  2d  off-diagonal  blocks  in  the 
case  of  simplex  elements.  These  blocks  represent  the  edges  of  the  neighboring  elements  of  one 
edge.  Each  block  is  dense  and  has  O  {m^  .p2(rf-i)^  entries.  Please  recall  that  p  is  the  polynomial 
degree  of  the  ansatz  functions,  d  is  the  spatial  dimension  of  the  domain  0,  and  m  is  the  number  of 
partial  differential  equations  (m  =  4  for  the  2-dimensional  Euler  or  Navier-Stokes  equations).  This 
structure  is  very  similar  to  that  of  a  normal  DG  discretization,  whereas  the  blocks  in  the  latter 
have  O  {m^  •  entries  and  thus  are  considerably  bigger  for  higher  polynomial  order  p. 

4  Adaptation  Procedure 

In  the  context  of  adjoint-based  (also  referred  to  as  target-  or  output-based)  error  estimation,  one 
is  interested  in  quantifying  the  error  of  a  specific  target  functional  Jh  '■  IP;  i-C- 

eh  ■=  Jh  {^)  -  Jh  i^h) ,  (46) 

where  x/j  is  the  approximation  to  x  in  X/^.  This  target  functional  can,  for  example,  represent  lift  or 
drag  coefficients  in  aerospace  applications.  Eor  the  derivation  of  the  adjoint-based  error  estimate 
we  expand  the  target  functional  in  a  Taylor  series  as  follows 

Jh  (5f)  -  Jh  i'^h)  =  J'h  [^/i]  (^  -^h)  +  0  (||x  -  .  (47) 

We  proceed  in  a  similar  manner  with  the  error  in  the  residual,  i.e. 

JJh  {^]yh)  -Afh  {^h]Jh)  =Afh  [^/i]  i^-^h;jh)  +  O  (||x  .  (48) 

As  our  discretization  is  consistent  the  first  term  Mh  (^^jh)  vanishes. 

Substituting  Eq.  (48)  into  Eq.  (47)  and  neglecting  the  quadratic  terms  yields 

Bh^rj  :=  -Afh  i^h]  ®/i)  (49) 

where  k/j  is  defined  by  the  so-called  adjoint  equation 

^h  (jh;  ^h)  =  J'h  (jh)  Jjh  e  Xft.  (50) 

The  adjoint  solution  'Zh  =  (^qh,yJh,  ^h^  £  ^h  represents  the  link  between  variations  in  the 
residual  and  in  the  target  functional. 

The  global  error  estimate  r]  can  then  be  restricted  to  a  single  element  to  yield  a  local  indicator 
to  drive  an  adaptation  procedure,  i.e. 

VK  ■=  |A4  {cih,Wh,Xh;'^h,Wh,0)  1^1 .  (51) 
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Here,  we  want  to  emphasize  that,  in  contrast  to  the  global  error  estimate,  we  ignore  the  contribution 
of  the  hybrid  adjoint  variable  A/^.  We  found  that  by  taking  the  whole  adjoint  into  account,  jumps 
across  element  interfaces  were  overly  penalized.  This  deserves  a  more  in-depth  analysis. 

Please  note,  that  the  functionals  Mh  and  and  their  jacobians  have  to  be  evaluated  in  a 
somewhat  richer  space  than  Xh,  namely  D  X/j.  Otherwise  the  weighted  residual  J\fh  {'^h','^h) 
would  be  identical  zero  as 

A/'/i  =  0  Vy/j  G  X/i.  (52) 

This  can  be  achieved  by  either  mesh  refinement  or  a  higher  polynomial  degree  of  the  ansatz  func¬ 
tions.  In  our  setting,  especially  when  using  a  hierarchical  basis,  the  latter  is  advantageous  with 
respect  to  implementational  effort  and  efficiency. 


4.1  Hybridization 

In  matrix  form,  the  adjoint  system  (see  Eq.  (50))  reads  as  follows 


■  A 

B 

R  ■ 

T 

Q 

’  F  ’ 

C 

D 

w 

= 

G 

L 

M 

N 

A 

H 

(53) 


Please  note,  that  in  our  formulation  H  =  0  as  A/i  is  not  defined  on  the  boundary  and  thus  the 
target  functional  depends  only  on  Wh  and  qh- 

As  the  overall  structure  of  the  adjoint  equation  is  similar  to  the  primal  system  (see  Eq.  (45)), 
one  can  also  apply  static  condensation  to  the  adjoint  system  which  then  yields  its  hybridized  form: 


-[  L  M 


■  A 

B  ' 

-1 

'  R  ' 

C 

D 

S 

A 


■  A 

B  ' 

-T 

'  F  ' 

C 

D 

G 

It  is  interesting  to  note  that  the  hybridized  adjoint  system  matrix  is  also  the  transpose  of  the 
hybridized  primal  system  matrix  (for  a  higher  polynomial  order).  This  is  very  beneficial  for  the 
implementation  as  the  routines  for  the  assembly  of  this  matrix  are  already  available. 

The  adjoint  solution  within  each  element  can  then  be  computed  with  the  aid  of  the  adjoint  local 
system,  given  by 


■  A 

B  ' 

T 

’  Q  ' 

F 

C 

D 

w 

G 

■  L  M  ]'^A 


(54) 


where  the  matrix  is  also  the  transpose  of  the  primal  local  matrix  (see  Eq.  (43)). 


4.2  Marking  Elements  for  Refinement 

After  having  obtained  a  localized  error  estimate,  we  choose  a  set  of  elements  to  be  refined.  The 
aim  of  our  marking  strategy  is  to  find  the  smallest  set  AI  C  7/i  such  that  the  error  contributed  by 
this  set  represents  a  certain  fraction  of  the  total  error,  i.e. 

(55) 

The  user-defined  parameter  6  is  of  course  problem  dependent.  It  can,  however,  be  tuned  for  a  big 
range  of  test  cases.  Please  note,  that  we  define  the  error  of  any  subset  of  Th  as  rjj^  :=  YIk^m  ^k- 
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4.3  Choosing  between  h-  and  p- Adaptation 


The  final  step  in  the  adaptation  procedure  is  the  decision  between  mesh  refinement  and  order 
enrichment.  On  each  element  a  smoothness  sensor  is  defined  as 

{w  —  W*,  W  —  w*)j^ 


Sk  ■= 


{w,  w] 


(56) 


K 


where  w*  is  the  element-wise  projection  of  w  to  the  next  smaller  polynomial  space,  given  by 


^  (57) 

Hence,  w  —  w*  represents  the  higher  order  components  of  the  solution  (see  Fig.  1).  As  we  use  an 
hierarchical  basis,  this  projection  is  very  cheap.  By  introducing  a  threshold  es,  a  decision  between 
mesh-refinement  and  p-enrichment  can  be  made,  i.e. 


Sk 


p-enrichment 

mesh-refinement 


(58) 


5  Results 

In  the  following  we  compare  our  in-house  HDG  and  DG  solvers  in  terms  of  degrees  of  freedom  and 
runtime.  The  DG  discretization  is  based  on  the  Lax-Friedrich  and  the  BR2  fluxes  for  convective  and 
viscous  terms,  respectively.  Both  solvers  share  the  same  computational  framework,  so  we  believe 
that  our  comparison  is  meaningful. 

We  apply  both  solvers  to  compressible  flow  problems,  including  inviscid  subsonic  and  transonic, 
and  subsonic  laminar  flow.  In  all  cases  we  show  results  for  pure  mesh-adaptation  (p  =  1 . . .  4)  and 
/ip-adaptation  (p  =  2  . . .  5). 

Please  note,  that  we  choose  the  same  relaxation  parameters  for  all  test  cases.  If  artificial 
viscosity  is  necessary,  we  use  eo  =  0.2  and  f3  =  0.  The  parameters  for  the  adaptation  procedure 
are  chosen  as  9  =  0.05  and  es  =  10“®.  This  set  of  parameters  seems  to  be  robust  for  both  DG  and 
HDG  for  a  broad  range  of  test  cases.  In  order  to  approximate  the  error  in  drag,  reference  solutions 
on  /i-p-adapted  meshes  with  more  than  2  •  10®  degrees  of  freedom  (please  note,  that  we  refer  to 
ndofu,  whenever  we  speak  of  degrees  of  freedom  in  the  following  as  this  is  a  good  measure  for  the 
resolution)  are  used. 

Before  we  turn  our  attention  to  the  adaptive  computations,  we  want  to  compare  runtimes  for 
both  methods  on  a  fixed  mesh  for  several  polynomial  orders.  This  way,  we  can  a  priori  learn  which 
improvement  can  be  expected.  In  Fig.  2  timings  for  the  assembly  procedure  and  the  iterative  solver 
are  given.  We  show  Euler  and  a  Navier-Stokes  computation  on  a  mesh  with  2396  elements  and  3544 
interior  faces.  We  used  polynomial  orders  from  p  =  0  to  p  =  6.  We  want  to  emphasize  that  the 
necessary  Newton  and  GMRES  iterations  for  both  HDG  and  DG  were  comparable.  As  there  are 
more  faces  than  elements,  DG  is  faster  than  HDG  for  p  =  0  and  p  =  1.  However,  already  for  p  =  2 
HDG  catches  up.  At  p  =  6  there  is  a  ratio  of  2.5  for  the  Euler  test  case  and  2.1  for  the  Navier-Stokes 
test  case.  For  the  Euler  test  case  it  is  interesting  to  note,  that  the  iterative  solver  dominates  the 
computational  time  for  DG.  Eor  HDG  it  is  the  other  way  around.  This  has  two  reasons:  firstly,  the 
assembly  is  more  involved  due  to  the  local  solves;  and  secondly,  the  global  system  is  considerably 
smaller  for  higher  polynomial  degree.  In  the  case  of  a  Navier-Stokes  computation,  the  DG  assembly 
takes  over  the  dominating  part  as  the  lifting  operators  are  very  expensive  to  compute.  The  time 
for  the  HDG  assembly  is  also  increased  which  is  among  other  things  due  to  the  introduction  of  the 
gradient.  For  both  HDG  and  DG,  the  time  spent  in  the  iterative  solver  is  comparable  for  both 
Euler  and  Navier-Stokes. 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


(a)  Coarse  mesh 


(b)  Adapted  mesh 

Figure  1:  Smoothness  sensor  for  a  transonic  test  case 
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(a)  Euler 


(b)  Navier-Stokes 


Figure  2:  Runtime  comparison  of  the  hybridized  and  non-hybridized  DG  method  for  a  fixed  mesh 
and  varying  polynomial  degree  (Solid:  HDG,  Shaded:  DG). 


5.1  Subsonic  Inviscid  Flow  over  the  NACA  0012  Airfoil 

In  the  first  test  case,  we  consider  subsonic  inviscid  flow  over  the  NACA  0012  airfoil  which  is  defined 
by 

y  =  ±0.6  (0.2969^/x  -  0.1260x  -  O.SSlGx^  ±  0.2843x^  -  O.lOlSx"^)  (59) 

with  X  G  [0, 1].  Using  this  definition,  the  airfoil  would  have  a  finite  trailing  edge  thickness  of  .252  %. 

In  order  to  obtain  a  sharp  trailing  edge  we  modify  the  coefficient,  i.e. 

y  =  ±0.6  (0.2969\/x  -  0.1260x  -  0.3516x2  ±  0.2843x^  -  0.1036x'‘)  .  (60) 

The  flow  is  characterized  by  a  free  stream  Mach  number  of  Maoo  =  0.5  and  an  angle  of  attack 

of  a  =  2°.  In  Fig.  3  the  baseline  mesh  for  the  Euler  test  cases  (subsonic  and  transonic)  can  be  seen. 
It  consists  of  719  elements  and  its  far  field  is  over  a  1000  chords  away. 

Admissible  target  functionals  defined  on  the  boundary  for  the  Euler  equations  are  given  by  the 
weighted  pressure  along  wall  boundaries,  i.e. 

J(w)  =  xjj  ■  (pn)  dcj  (61) 

Jdn 

where  n  is  the  outward  pointing  normal.  By  using  xf)  =  (cos  a,  sin  a)^  or  ■?/>  =  (—  sin  a,  cos  a)^ 

along  wall  boundaries  and  0  otherwise,  the  functional  represents  the  pressure  drag  coefficient  cd 
or  the  pressure  lift  coefficient  cl,  respectively.  Coo  is  a  normalized  reference  value  defined  by 
Coo  =  ^ylVIa^poo^-  Here,  I  is  the  chord  length  of  the  airfoil. 

In  Fig.  4a,  a  purely  /i-adapted  mesh  can  be  seen.  The  most  refined  regions  are  the  leading 
and  trailing  edge.  The  first  is  of  importance  as  the  flow  experiences  high  gradients  towards  the 
stagnation  point.  Refinement  of  the  latter  is  necessary  due  to  the  sharp  trailing  edge  and  the 
slip-wall  boundary  conditions.  As  soon  as  the  error  in  these  two  regions  is  sufficiently  low,  other 
elements  close  to  the  airfoil  get  refined  as  well.  For  the  /ip-adapted  mesh  (see  Fig.  4b)  the  leading 
and  trailing  edge  are  refined  as  well.  All  other  regions,  however,  undergo  mostly  p-enrichment. 
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Figure  3;  Baseline  mesh  with  719  elements  for  inviscid  computations 


p 

1 

2 

3 

4 

hp 

^dg/^hdg 

1.293 

3.799 

3.551 

3.092 

2.209 

?^nz,DG/’^nz,HDG 

1.123 

2.150 

3.248 

4.281 

4.601 

Table  1:  Runtime  and  nonzero  ratios  for  a  fixed  error  level  (Maoo  =  0.5,  a  =  2°) 


In  terms  of  degrees  of  freedom,  both  HDG  and  DG  show  similar  results.  For  all  computations 
it  takes  some  adaptations  until  the  critical  regions,  leading  and  trailing  edge,  are  resolved.  From 
this  point  on,  one  can  see  the  benefit  of  a  higher  order  discretization:  the  error  drops  significantly 
faster  with  respect  to  degrees  of  freedom  and  computational  time  (see  Fig.  5  and  6). 

In  Tbl.  1,  we  give  the  runtime  ratios  for  a  fixed  error  level  (we  always  choose  the  minimum 
level  attained).  For  p  =  1,  HDG  and  DG  are  comparable  in  both  runtime  and  nonzero  entries. 
From  p  =  2  on,  HDG  is  faster  than  DG  and  needs  less  memories  for  the  global  system.  Here,  it 
is  important  to  note  that  the  adjoint  is  approximated  with  p  +  1,  so  that  already  for  lower  orders 
HDG  is  faster. 

5.2  Transonic  Inviscid  Flow  over  the  NACA  0012  Airfoil 

Next,  we  turn  our  attention  to  transonic  flow  which  develops  more  complex  features  (e.g.  compres¬ 
sion  shocks)  compared  to  the  subsonic  regime.  The  flow  is  characterized  by  a  free  stream  Mach 
number  of  Maoo  =  0.8  and  an  angle  of  attack  of  a  =  1.25°. 

In  Fig.  7a  a  purely  /i-adapted  mesh  can  be  seen.  The  adjoint  sensor  detects  all  regions  of 
relevance  for  the  drag:  the  upper  shock,  the  leading  and  trailing  edge,  and  the  lower  weak  shock. 
Further  refinement  is  added  upstream  of  the  shock,  where  the  adjoint  has  steep  gradients  and  thus 
needs  higher  resolution.  In  the  case  of  /ip-adaptation,  the  mesh-refinement  is  stronger  confined  to 
the  shock  region  and  the  trailing  edge.  The  other  features  undergo  p-enrichment. 

As  expected,  both  methods  show  a  similar  accuracy  for  a  given  number  of  degrees  of  free¬ 
dom.  The  computations  with  p  =  2 ...  4  outperform  p  =  1  but  are  comparable  to  each  other, 
/ip-adaptation  shows  very  good  results  which  is  due  to  the  accurate  prediction  of  the  solution 
smoothness  (see  Fig.  8  and  9). 

For  this  test  case,  HDG  is  faster  than  DG  from  p  =  1  on  (see  Tbl.  2).  The  hp-adaptive  run 
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(a)  Pure  /i-adapation  {p  =  2) 


(b)  /ip-adapation  (p  =  2  . . .  5) 

Figure  4:  Adapted  meshes  for  the  subsonic  Euler  test  case  (Maoo  =  0.5,  a  =  2°) 
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(a)  HDG  (b)  DG 

Figure  5:  Drag  convergence  with  respect  to  degrees  of  freedom  (Maoo  =  0.5,  a  =  2°) 
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3 

4 

hp 

tDc/tuDG 

2.636 

1.607 

2.154 

2.362 

2.628 

?T'nz,DG/’^nz,HDG 

1.237 

1.807 

3.098 

4.437 

2.411 

Table  2:  Runtime  and  nonzero  ratios  for  a  fixed  error  level  (Maoo  =  0.8,  a  =  1.25°) 


is  more  than  2.5  times  faster.  The  ratio  of  necessary  nonzero  entries  attains  its  highest  value  for 
p  =  A.  In  the  case  of  /ip-adaptation,  this  ratio  is  not  as  high,  as  in  the  shock  region  a  lot  of  elements 
with  p  =  2  exist. 

5.3  Subsonic  Laminar  Flow  over  the  NACA  0012  Airfoil 

Finally,  we  consider  viscous  flow  in  the  subsonic  regime.  The  free  stream  Mach  number  is  Maoo  = 
0.5,  the  angle  of  attack  a  =  1°  and  the  Reynolds  number  Re  =  5000.  Due  to  the  latter,  a  thin 
boundary  layer  develops  around  the  airfoil. 

The  baseline  mesh  for  the  Navier-Stokes  test  case  is  more  refined  around  the  airfoil  such  that 
the  boundary  layer  is  correctly  captured  (see  Fig.  10).  It  consists  of  1781  elements  and  its  far  field 
is  over  a  1000  chords  away. 

Admissible  target  functionals  defined  on  the  boundary  for  the  Navier-Stokes  equations  are  given 
by  the  weighted  boundary  flux  along  wall  boundaries,  i.e. 

J{w,'Vw)=  /  '0  •  (pn  —  rn)  d(T  (62) 

Jdn 

where  n  is  the  outward  pointing  normal.  Here,  0  is  nonzero  only  on  wall  boundaries.  By  using 
0  =  ^  (cos  a,  sin  a)^  or  0  =  ^  (— sin  a,  cos  a)^  along  wall  boundaries  and  otherwise  0,  the 
functional  represents  the  viscous  drag  coefficient  cd  or  the  viscous  lift  coefficient  cl,  respectively. 

Both  the  h-adapted  mesh  (see  Fig.  11a)  and  the  /ip-adapted  mesh  (see  Fig.  11b)  undergo 
refinement  within  the  boundary  layer  and  the  wake  region.  The  mesh  refinement  for  the  hp- 
adaptive  run  is  however  more  restricted  to  the  leading  edge  region  where  the  boundary  layer 
develops.  Further  downstream,  p-enrichment  is  used  as  soon  as  the  necessary  mesh-resolution 
is  reached. 
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(a)  HDG  (b)  DG 

Figure  6:  Drag  convergence  with  respect  to  time  (Maoo  =  0.5,  a  =  2°) 


p 

1 

2 

3 

4 

hp 

iDG/^HDG 

2.244 

2.572 

2.397 

2.288 

3.122 

?^nz,DG/?^•nz,HDG 

1.196 

2.235 

3.472 

4.531 

3.463 

Table  3:  Runtime  and  nonzero  ratios  for  a  fixed  error  level  (Maoo  =  0.5,  a  =  1°,  Re  =  5000) 


In  terms  of  accuracy  versus  degrees  of  freedom,  HDG  and  DG  perform  comparably  well.  The 
higher  the  polynomial  degree  the  more  accurate  and  efficient  the  computations  are  for  both  HDG 
and  DG  (see  Fig.  12  and  13).  The  difference  between  /ip,  p  =  3  and  p  =  4  is  not  as  big,  though. 
This  might  lead  to  the  conclusion  that  isotropic  mesh  refinement  is  not  longer  efficiently  applicable 
in  cases  involving  strong  gradients.  Hence,  the  efficiency  of  the  adaptation  procedure  is  rather 
limited  by  the  mesh  refinement  strategy. 

Concerning  the  timings,  we  can  see  a  similar  trend  as  in  the  previous  test  cases  (see  Tbl.  3). 
For  p  =  1 . . .  4,  HDG  is  more  than  twice  as  fast.  The  /ip-adaptive  HDG  computation  is  even  three 
times  as  fast  compared  to  the  DG  run.  From  p  =  2  on,  the  savings  in  nonzero  entries  for  HDG 
become  significant. 
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(b)  /ip-adaptation  (p  =  2  . . .  5) 

Figure  7;  Adapted  meshes  for  the  transonic  Euler  test  case  (Maoo  =  0.8,  a  =  1.25°) 
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(a)  HDG  (b)  DG 

Figure  8:  Drag  convergence  with  respect  to  degrees  of  freedom  (Maoo  =  0.8,  a  =  1.25°) 


t/to 

(a)  HDG 
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(b)  DG 


Figure  9:  Drag  convergence  with  respect  to  time  (Maoo  =  0.8,  a  =  1.25°) 
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Figure  10:  Baseline  mesh  with  1781  elements  for  viscous  computations 
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(a)  Pure  ^-adaptation  {p  =  2) 


(b)  /ip-adaptation  (p  =  2  . . .  5) 

Figure  11;  Adapted  meshes  for  the  Navier-Stokes  test  case  (Maoo  =  0.5,  a  =  1°,  Re  =  5000) 
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(a)  HDG  (b)  DG 


5ure  12:  Drag  convergence  with  respect  to  degrees  of  freedom  (Maoo  =  0.5,  a  =  1°,  Re  =  5000) 


t/to  t/to 

(a)  HDG  (b)  DG 

Figure  13:  Drag  convergence  with  respect  to  time  (Maoo  =  0.5,  a  =  1°,  Re  =  5000) 
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